# Replication code for Di Landro, Gonzalo (2024) Party Behavior and the Gender Voting Gap
# British Journal of Political Science

library(ggplot2)
library(stargazer)
library(dplyr)
library(readxl)
library(lme4)
library(interplot)
library(gridExtra)
library(mlogit)


# Loading data ------------------------------------------------------------
load("gap_dataset.RData")
gap_left<- subset(gap,gap$parfam<31)

#Keeping most voted per family
gap_left <- gap_left %>%
  group_by(country, date, parfam) %>%
  filter(pervote == max(pervote, na.rm = TRUE)) %>%
  ungroup()  

# Figure 1 ----------------------------------------------------------------------

# A
A <- ggplot(gap_left, aes(as.numeric(as.character(eyear)), F_M)) +
  geom_smooth(method = "loess", se = FALSE, color = "#F8766D") +
  theme_bw() +
  labs(
    title = "A. Pooled sample",
    y = "F/M party voter ratio",
    x = ""
  ) +
  geom_point( color = "#F8766D",size=0.2)+

  theme(
    legend.position = "none",
    plot.title = element_text(size = 12, face = "bold"),
    axis.text.x = element_text(angle = 40, hjust = 1, size = 8.5),
    strip.text = element_text(size = 12)
  ) +  
  geom_hline(yintercept = 1, linetype = "dashed") +
  geom_smooth(
    aes(y = as.numeric(flfp) / 25),
    colour = "black",
    size = 1.5,
    method = "loess",
    linetype = "dotted",
    se = FALSE
  ) +
  scale_y_continuous(
    limits = c(0, 3.5),
    breaks = seq(0, 3.5, by = 0.5),
    sec.axis = sec_axis(~ . * 25, name = "", breaks = seq(0,100, by = 20))
  ) +
  scale_x_continuous(
    limits = c(1985, 2020),
    breaks = seq(1980, 2020, by = 10)
  ) +
  theme(
    axis.text = element_text(size = 10),
    axis.title = element_text(size = 13),
    plot.title = element_text(size = 16, face = "bold")
  )

# B
gap_left$country_toplot <- gap_left$country
gap_left[which(gap_left$country == "United Kingdom"), "country_toplot"] <- "UK"

B <- ggplot(gap_left, aes(as.numeric(as.character(eyear)), F_M)) +
  geom_smooth(method = "loess", se = FALSE, color = "#F8766D") +
  theme_bw() +
  labs(
    title = "B. Country samples",
    y = "",
    x = ""
  ) +
  theme(
    legend.position = "none",
    plot.title = element_text(size = 12, face = "bold"),
    axis.text.x = element_text(angle = 40, hjust = 1, size = 8.5),
    strip.text = element_text(size = 12)
  ) +  
  facet_wrap(~country_toplot) +
  geom_hline(yintercept = 1, linetype = "dashed") +
  geom_line(aes(y = as.numeric(flfp) / 25), colour = "black") +
  geom_point(
    aes(y = as.numeric(flfp) / 25),
    colour = "black",
    size = 0.5,
    shape = 22,
    fill = "black"
  ) +
  scale_y_continuous(
    sec.axis = sec_axis(~ . * 25, name = "Female Labor Force Participation rate (FLFP)",breaks = seq(0,100, by = 20))
  ) +
  theme(
    axis.text = element_text(size = 10),
    axis.title = element_text(size = 13),
    plot.title = element_text(size = 16, face = "bold")
  ) +
  scale_x_continuous(
    limits = c(1985, 2020),
    breaks = seq(1980, 2020, by = 10)
  )

# Figure
grid.arrange(A, B, ncol = 2)


# Table 1 -----------------------------------------------------------------

# Model m0
m0 <- lmer(log(F_M) ~ flfp_lag1 + divorce_rate_lag1 + per_chr + (1 | party) + (1 | country_year), 
           data = gap_left, REML = FALSE)

# Model m2
m2 <- lmer(log(F_M) ~ flfp_lag1 + gender_labor_equality + divorce_rate_lag1 + per_chr + share_women_lag1 + 
             female_leader_lag1 + (1 | party) + (1 | country_year), 
           data = gap_left, REML = FALSE)

# Model m3
m3 <- lmer(log(F_M) ~ flfp_lag1 * gender_labor_equality + divorce_rate_lag1 + per_chr + share_women_lag1 + 
             female_leader_lag1 + (1 | party) + (1 | country_year), 
           data = gap_left, REML = FALSE)

# Variable order for stargazer
order <- c("flfp_lag1", "gender_labor_equality", "share_women_lag1", "female_leader_lag1", 
           "divorce_rate_lag1", "per_chr")

# Random effects for m0
RE0 <- as.data.frame(VarCorr(m0))[, 4]
RE0 <- round(RE0, 3)

# Random effects for m2
RE2 <- as.data.frame(VarCorr(m2))[, 4]
RE2 <- round(RE2, 3)

# Random effects for m3
RE3 <- as.data.frame(VarCorr(m3))[, 4]
RE3 <- round(RE3, 3)

# Table
stargazer(m0, m2, m3,
          title = "Table 1. Determinants of the gender voting gap",
          add.lines = list(
            c("Var country-year", RE0[1], RE2[1], RE3[1]),
            c("Var party", RE0[2], RE2[2], RE3[2]),
            c("Var residual", RE0[3], RE2[3], RE3[3]),
            c("N country-years", nrow(ranef(m0)[[1]]), nrow(ranef(m2)[[1]]), nrow(ranef(m3)[[1]])),
            c("N parties", nrow(ranef(m0)[[2]]), nrow(ranef(m2)[[2]]), nrow(ranef(m3)[[2]]))
          ),
          order = order,
          star.cutoffs = c(0.05, 0.01, 0.001),
          dep.var.labels = 'log(F/M) party voter ratio',
          single.row = TRUE,
          type = "text")

# Figure 2 ----------------------------------------------------------------

interplot(m = m3, var1 = "flfp_lag1", var2 = "gender_labor_equality", hist = TRUE, ci = 0.95) + 
  theme_bw() +
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(
    title = "",
    y = "Effect of FLFP on log(F/M)",
    x = "Parties' support for gender labor equality"
  ) +
  theme(
    axis.text = element_text(size = 10),
    axis.title = element_text(size = 11),
    plot.title = element_text(size = 16, face = "bold"),
    plot.caption = element_text(hjust = 0, size = 8)
  )

# Figure A1: Parties’ support for gender labor equality  --------------


# Load the dataset into gle_averages
gle_averages <- gap

# Filter the dataset to keep only the party with the highest percentage of votes (pervote) for each country, date, and party family
gle_averages <- gle_averages %>%
  group_by(country, date, parfam) %>%
  filter(pervote == max(pervote, na.rm = TRUE)) %>%
  ungroup()  

# Calculate the mean gender labor equality (GLE) score for each party family
gle_averages <- gle_averages %>%
  group_by(parfam) %>%
  summarize(mean_gle = mean(gender_labor_equality, na.rm = TRUE))

# Create a factor variable indicating left-wing parties (parfam < 31)
gle_averages$left_parties <- as.factor(ifelse(gle_averages$parfam < 31, 1, 0))

# Set the seed for reproducibility
set.seed(1)

# Plot A: Bar plot of average GLE scores by party family
A2a <- ggplot(gle_averages, aes(x = factor(parfam), y = mean_gle, fill = left_parties)) +
  geom_bar(stat = "identity") +
  theme_bw() +
  labs(
    title = "A. Parties' support for gender labor equality by family",
    y = "Average score",
    x = ""
  ) +
  scale_fill_manual(values = c("#619CFF", "#F8766D")) +
  theme(
    axis.text = element_text(size = 10),
    legend.position = "none",
    axis.title = element_text(size = 10),
    plot.title = element_text(size = 10, face = "bold"),
    plot.caption = element_text(hjust = 1, size = 8)
  ) +
  scale_x_discrete(
    labels = c(
      "Ecological", "Socialist", "Social Democratic", "Liberal",
      "Christian Democratic", "Conservative", "Nationalist", "Agrarian", "Ethnic", "Special-Issue"
    )
  ) +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
    legend.position = "none"
  )

# Plot B: Jitter plot of GLE scores by country for left-wing parties
A2b <- ggplot(gap_left, aes(x = factor(country), y = gender_labor_equality, colour = as.factor(parfam))) +
  geom_jitter(width = 0.25, size = 0.5) +
  theme_bw() +
  labs(
    title = "B. Left-wing parties' support for gender labor equality by country",
    y = "Score",
    x = "",
    caption = ""
  ) +
  theme(
    axis.text = element_text(size = 10),
    legend.title = element_blank(),
    legend.position = "right",
    axis.title = element_text(size = 10),
    plot.title = element_text(size = 10, face = "bold"),
    plot.caption = element_text(hjust = 1, size = 8),
    axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)
  ) +
  scale_color_manual(values = c("#008F2A", "#8049CC", "#BF1B0B"),
                     labels = c('Ecological', 'Socialist', 'Social Democratic'))

# Figure
grid.arrange(A2a, A2b, ncol = 1, nrow = 2)


# Figure A2:  Left-wing parties’ support for gender labor equality in Western Europe by country ---------------------------------------------------------------
ggplot(gap_left, aes(as.numeric(as.character(eyear)), gender_labor_equality)) +
  geom_smooth(method = "loess", se = FALSE) +
  facet_wrap(~ country) +
  scale_color_manual(values = c('Black')) +
  theme_bw() +
  labs(
    title = "",
    y = "Score",
    x = "",
    caption = ""
  ) +
  theme(
    axis.text = element_text(size = 10),
    legend.position = "none",
    axis.title = element_text(size = 10),
    axis.text.x = element_text(angle = 40, hjust = 1),
    plot.title = element_text(size = 10, face = "bold"),
    plot.caption = element_text(hjust = 1, size = 8)
  )

# Figure A3: Quartiles of left-wing parties’ support for gender labor equality in Western Europe ---------------------------------------------------------------

# mean and quartiles
mean_gle <- mean(gap_left$gender_labor_equality, na.rm = TRUE)
quartiles <- quantile(gap_left$gender_labor_equality, probs = c(0.25, 0.5, 0.75), na.rm = TRUE)

ggplot(gap_left) + 
  geom_boxplot(aes(y = gender_labor_equality, fill = "red")) + 
  scale_x_discrete() +
  theme_bw() +
  coord_flip() +
  labs(
    title = "",
    y = "Score",
    x = "",
    caption = ""
  ) +
  theme(
    axis.text = element_text(size = 10),
    legend.position = "none",
    axis.title = element_text(size = 10),
    plot.title = element_text(size = 10, face = "bold"),
    plot.caption = element_text(hjust = 1, size = 8)
  )

# Figure A4: Correlation between CMP and VParty indicators ---------------------------------------------------------------

correlation <- cor(gap_left$per504, gap_left$gender_labor_equality, use = "pairwise.complete.obs") 
print(correlation)  # corr = 0.15

ggplot(gap_left, aes(x = per504, y = gender_labor_equality)) +
  geom_point() +
  geom_smooth(method = "lm") +
  labs(
    title = "",
    y = "Parties' support for gender labor equality (V-Party)",
    x = "Welfare expansion (Per504 CMP)",
    caption = ""
  ) +
  theme_bw() +
  theme(
    axis.text = element_text(size = 10),
    legend.position = "none",
    axis.title = element_text(size = 10),
    plot.title = element_text(size = 10, face = "bold"),
    plot.caption = element_text(hjust = 1, size = 8)
  ) +
  geom_label(aes(label = "Corr = 0.15", x = 31, y = 0.65), size = 4)


# Figure A5: Distribution of F/M and log (F/M) ---------------------------------------

# density of F/M
p1 <- ggplot(gap_left, aes(F_M)) +
  geom_density() +
  theme_bw() +
  labs(
    title = "A",
    y = "Density",
    x = "F/M",
    caption = ""
  ) +
  theme(
    axis.text = element_text(size = 10),
    legend.position = "none",
    axis.title = element_text(size = 10),
    plot.title = element_text(size = 10, face = "bold"),
    plot.caption = element_text(hjust = 1, size = 8)
  )

# density of log(F/M)
p2 <- ggplot(gap_left, aes(log(F_M))) +
  geom_density() +
  theme_bw() +
  labs(
    title = "B",
    y = "Density",
    x = "log(F/M)",
    caption = ""
  ) +
  theme(
    axis.text = element_text(size = 10),
    legend.position = "none",
    axis.title = element_text(size = 10),
    plot.title = element_text(size = 10, face = "bold"),
    plot.caption = element_text(hjust = 1, size = 8)
  )

# Figure 
grid.arrange(p1, p2, ncol = 2)

# Figure A6: Deciles of FLFP ---------------------------------------------------------------

# decile groups
deciles <- quantile(gap_left$flfp_lag1, probs = seq(0.1, 0.9, 0.1), na.rm = TRUE)
deciles_df <- stack(deciles)
names(deciles_df) <- c("value", "decile")
deciles_vec <- deciles_df$value

# Function to color
StatAreaUnderDensity <- ggproto(
  "StatAreaUnderDensity", Stat,
  required_aes = "x",
  compute_group = function(data, scales, xlim = NULL, n = 50) {
    fun <- approxfun(density(data$x))
    StatFunction$compute_group(data, scales, fun = fun, xlim = xlim, n = n)
  }
)

stat_aud <- function(mapping = NULL, data = NULL, geom = "area",
                     position = "identity", na.rm = FALSE, show.legend = NA, 
                     inherit.aes = TRUE, n = 50, xlim = NULL, ...) {
  layer(
    stat = StatAreaUnderDensity, data = data, mapping = mapping, geom = geom, 
    position = position, show.legend = show.legend, inherit.aes = inherit.aes,
    params = list(xlim = xlim, n = n, ...)
  )
}

# Figure
ggplot(gap_left, aes(flfp_lag1)) +
  geom_density() +
  theme_bw() +
  geom_vline(xintercept = deciles, linetype = "dashed", color = "red") +
  labs(
    title = "",
    y = "Density",
    x = "FLFP_lag",
    caption = ""
  ) +
  theme(
    axis.text = element_text(size = 10),
    legend.position = "none",
    axis.title = element_text(size = 10),
    plot.title = element_text(size = 10, face = "bold"),
    plot.caption = element_text(hjust = 1, size = 8)
  ) +
  stat_aud(
    geom = "area",
    aes(fill = flfp_lag1),
    xlim = c(33.489, 46.33333),
    alpha = 0.2
  )


# Table A1: Models with random effects for country--------------------------------------------


# Model m0
m0 <- lmer(log(F_M) ~ flfp_lag1 + divorce_rate_lag1 + per_chr + prop2 +
             (1 | party) + (1 | country) + (1 | country_year), 
           data = gap_left, REML = FALSE)

# Model m2
m2 <- lmer(log(F_M) ~ flfp_lag1 + gender_labor_equality + divorce_rate_lag1 + per_chr + 
             share_women_lag1 + female_leader_lag1 + prop2 + (1 | party) + (1 | country) + (1 | country_year), 
           data = gap_left, REML = FALSE)

# Model m3
m3 <- lmer(log(F_M) ~ flfp_lag1 * gender_labor_equality + divorce_rate_lag1 + per_chr + 
             share_women_lag1 + female_leader_lag1 + prop2 + (1 | party) + (1 | country) + (1 | country_year), 
           data = gap_left, REML = FALSE)

# Variable order for stargazer
order <- c("flfp_lag1", "gender_labor_equality", "share_women_lag1", "female_leader_lag1", 
           "divorce_rate_lag1", "per_chr", "prop2")

# Extract and round random effects for each model
RE0 <- round(as.data.frame(VarCorr(m0))[, 4], 3)
RE2 <- round(as.data.frame(VarCorr(m2))[, 4], 3)
RE3 <- round(as.data.frame(VarCorr(m3))[, 4], 3)

# Table
stargazer(m0, m2, m3,
          title = "Models with random effects for country",
          add.lines = list(
            c("Var country-year", RE0[1], RE2[1], RE3[1]),
            c("Var party", RE0[2], RE2[2], RE3[2]),
            c("Var country", RE0[3], RE2[3], RE3[3]),
            c("Var residual", RE0[4], RE2[4], RE3[4]),
            c("N country-years", nrow(ranef(m0)[[1]]), nrow(ranef(m2)[[1]]), nrow(ranef(m3)[[1]])),
            c("N parties", nrow(ranef(m0)[[2]]), nrow(ranef(m2)[[2]]), nrow(ranef(m3)[[2]])),
            c("N countries", nrow(ranef(m0)[[3]]), nrow(ranef(m2)[[3]]), nrow(ranef(m3)[[3]]))
          ),
          order = order,
          dep.var.labels = 'log F/M party voter ratio',
          single.row = TRUE,
          star.cutoffs = c(0.05, 0.01, 0.001),
          type = "text")


# Table A2: Models controlling for right-wing parties’ (RWPs) behavior -----------------

# Subset for right-wing parties
gap_right <- subset(gap, gap$parfam > 30 & gap$parfam < 71)

# Keeping right-wing party with highest gender labor equality
gap_right <- gap_right %>%
  group_by(country, date) %>%
  filter(gender_labor_equality == max(gender_labor_equality, na.rm = TRUE)) %>%
  ungroup()

# Select relevant columns and rename them
gap_right <- gap_right[, c("country", "date", "party", "share_women_lag1", "female_leader_lag1", "gender_labor_equality")]
colnames(gap_right)[3:6] <- c("partyR", "share_women_lag1R", "female_leader_lag1R", "gender_labor_equalityR")

# Merge with left-wing party data
gap_combined <- merge(gap_left, gap_right, by = c("country", "date"), all.x = TRUE)

# Models
m0 <- lmer(log(F_M) ~ flfp_lag1 + divorce_rate_lag1 + per_chr + 
             (1 | party) + (1 | country_year), 
           data = gap_combined, REML = FALSE)

m2 <- lmer(log(F_M) ~ flfp_lag1 + divorce_rate_lag1 + per_chr + 
             share_women_lag1 + female_leader_lag1 + gender_labor_equality + 
             share_women_lag1R + female_leader_lag1R + gender_labor_equalityR + 
             (1 | party) + (1 | country_year), 
           data = gap_combined, REML = FALSE)

m3 <- lmer(log(F_M) ~ flfp_lag1 * gender_labor_equality + divorce_rate_lag1 + per_chr + 
             share_women_lag1 + female_leader_lag1 + 
             share_women_lag1R + female_leader_lag1R + gender_labor_equalityR + 
             (1 | party) + (1 | country_year), 
           data = gap_combined, REML = FALSE)


# Variable order for stargazer
order <- c("flfp_lag1", "divorce_rate_lag1", "per_chr", "gender_labor_equality", 
           "share_women_lag1", "female_leader_lag1", "share_women_lag1R", "female_leader_lag1R", 
           "gender_labor_equalityR")


# Extract and round random effects for each model
RE0 <- round(as.data.frame(VarCorr(m0))[, 4], 3)
RE2 <- round(as.data.frame(VarCorr(m2))[, 4], 3)
RE3 <- round(as.data.frame(VarCorr(m3))[, 4], 3)

# Table
stargazer(m0, m2, m3,
          title = "Models controlling for non-left parties' behavior",
          add.lines = list(
            c("Var country-year", RE0[1], RE2[1], RE3[1]),
            c("Var party", RE0[2], RE2[2], RE3[2]),
            c("Var residual", RE0[3], RE2[3], RE3[3]),
            c("N country-years", nrow(ranef(m0)[[1]]), nrow(ranef(m2)[[1]]), nrow(ranef(m3)[[1]])),
            c("N parties", nrow(ranef(m0)[[2]]), nrow(ranef(m2)[[2]]), nrow(ranef(m3)[[2]]))
          ),
          order = order,
          star.cutoffs = c(0.05, 0.01, 0.001),
          dep.var.labels = 'log(F/M) party voter ratio',
          single.row = TRUE,
          type = "text")


# Table A3: Models controlling for the share of the rural population ------------------------------------------------------------

# Model m0
m0 <- lmer(log(F_M) ~ flfp_lag1 + divorce_rate_lag1 + per_chr + per_rural + 
             (1 | party) + (1 | country_year), 
           data = gap_left, REML = FALSE)

# Model m2
m2 <- lmer(log(F_M) ~ flfp_lag1 + gender_labor_equality + divorce_rate_lag1 + per_chr + 
             per_rural + share_women_lag1 + female_leader_lag1 + 
             (1 | party) + (1 | country_year), 
           data = gap_left, REML = FALSE)

# Model m3
m3 <- lmer(log(F_M) ~ flfp_lag1 * gender_labor_equality + divorce_rate_lag1 + per_chr + 
             per_rural + share_women_lag1 + female_leader_lag1 + 
             (1 | party) + (1 | country_year), 
           data = gap_left, REML = FALSE)

# Variable order for stargazer
order <- c("flfp_lag1", "gender_labor_equality", "share_women_lag1", "female_leader_lag1", 
           "divorce_rate_lag1", "per_chr", "per_rural")

# Extract and round random effects for each model
RE0 <- round(as.data.frame(VarCorr(m0))[, 4], 3)
RE2 <- round(as.data.frame(VarCorr(m2))[, 4], 3)
RE3 <- round(as.data.frame(VarCorr(m3))[, 4], 3)

# Table
stargazer(m0, m2, m3,
          title = "Models controlling for the share of rural population",
          add.lines = list(
            c("Var country-year", RE0[1], RE2[1], RE3[1]),
            c("Var party", RE0[2], RE2[2], RE3[2]),
            c("Var residual", RE0[3], RE2[3], RE3[3]),
            c("N country-years", nrow(ranef(m0)[[1]]), nrow(ranef(m2)[[1]]), nrow(ranef(m3)[[1]])),
            c("N parties", nrow(ranef(m0)[[2]]), nrow(ranef(m2)[[2]]), nrow(ranef(m3)[[2]]))
          ),
          order = order,
          star.cutoffs = c(0.05, 0.01, 0.001),
          dep.var.labels = 'log(F/M) party voter ratio',
          single.row = TRUE,
          type = "text")




# Table A4: Models excluding outliers --------------------------------------------------------

# Subset for left-wing parties
gap_outliers <- subset(gap, gap$parfam < 31)

# Keep only the party with the highest percentage of votes (pervote) for each country, date, and party family
gap_outliers <- gap_outliers %>%
  group_by(country, date, parfam) %>%
  filter(pervote == max(pervote, na.rm = TRUE)) %>%
  ungroup()


# mean F/M
mean_FM <- mean(gap_outliers$F_M, na.rm = TRUE)

# Calculate the standard deviation
sd_FM <- sd(gap_outliers$F_M, na.rm = TRUE)

# Calculate two standard deviations above and below the mean
two_sd_above_FM <- mean_FM + 2 * sd_FM
two_sd_below_FM <- mean_FM - 2 * sd_FM

# Remove specific outliers and save as new dataset
gap_outliers_v2 <- subset(gap_outliers, F_M > two_sd_below_FM & 
                            F_M < two_sd_above_FM)

# Remove gender labor equality outliers and F_M
gap_outliers_v2 <- subset(gap_outliers_v2, !(country_year == "Italy 1994" & party == 32220)) # DPL 
gap_outliers_v2 <- subset(gap_outliers_v2, !(country_year == "Italy 1996" & party == 32220)) # DPL 
gap_outliers_v2 <- subset(gap_outliers_v2, !(country_year == "Italy 2006" & party == 32110)) # GF 


# Remove  only  gender labor equality outliers
gap_outliers <- subset(gap_outliers, !(country_year == "Italy 1994" & party == 32220)) # DPL 
gap_outliers <- subset(gap_outliers, !(country_year == "Italy 1996" & party == 32220)) # DPL 
gap_outliers <- subset(gap_outliers, !(country_year == "Italy 2006" & party == 32110)) # GF 


# Model m0
m0 <- lmer(log(F_M) ~ flfp_lag1 + divorce_rate_lag1 + per_chr + 
             (1 | party) + (1 | country_year), 
           data = gap_outliers, REML = FALSE)

# Model m2
m2 <- lmer(log(F_M) ~ flfp_lag1 + gender_labor_equality + divorce_rate_lag1 + per_chr + 
             share_women_lag1 + female_leader_lag1 + 
             (1 | party) + (1 | country_year), 
           data = gap_outliers, REML = FALSE)

# Model m3
m3 <- lmer(log(F_M) ~ flfp_lag1 * gender_labor_equality + divorce_rate_lag1 + per_chr + 
             share_women_lag1 + female_leader_lag1 + 
             (1 | party) + (1 | country_year), 
           data = gap_outliers, REML = FALSE)

# Model m4
m4 <- lmer(log(F_M) ~ flfp_lag1 * gender_labor_equality + divorce_rate_lag1 + per_chr + 
             share_women_lag1 + female_leader_lag1 + 
             (1 | party) + (1 | country_year), 
           data = gap_outliers_v2, REML = FALSE)

# Variable order for stargazer
order <- c("flfp_lag1", "gender_labor_equality", "share_women_lag1", "female_leader_lag1", 
           "divorce_rate_lag1", "per_chr")

# Extract and round random effects for each model
RE0 <- round(as.data.frame(VarCorr(m0))[, 4], 3)
RE2 <- round(as.data.frame(VarCorr(m2))[, 4], 3)
RE3 <- round(as.data.frame(VarCorr(m3))[, 4], 3)
RE4 <- round(as.data.frame(VarCorr(m4))[, 4], 3)

# Table
stargazer(m0, m2, m3, m4,
          title = "Excluding outliers",
          add.lines = list(
            c("Var country-year", RE0[1], RE2[1], RE3[1], RE4[1]),
            c("Var party", RE0[2], RE2[2], RE3[2], RE4[2]),
            c("Var residual", RE0[3], RE2[3], RE3[3], RE4[3]),
            c("N country-years", nrow(ranef(m0)[[1]]), nrow(ranef(m2)[[1]]), nrow(ranef(m3)[[1]]), nrow(ranef(m4)[[1]])),
            c("N parties", nrow(ranef(m0)[[2]]), nrow(ranef(m2)[[2]]), nrow(ranef(m3)[[2]]), nrow(ranef(m4)[[2]]))
          ),
          order = order,
          star.cutoffs = c(0.05, 0.01, 0.001),
          dep.var.labels = 'log(F/M) party voter ratio',
          single.row = TRUE,
          type = "text")

# Interaction plot
interplot(m = m4, var1 = "flfp_lag1", var2 = "gender_labor_equality", hist = TRUE, ci = 0.95) + 
  theme_bw() +
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(
    title = "",
    y = "Effect of FLFP on log(F/M)",
    x = "Parties' support for gender labor equality"
  ) +
  theme(
    axis.text = element_text(size = 10),
    axis.title = element_text(size = 11),
    plot.title = element_text(size = 16, face = "bold"),
    plot.caption = element_text(hjust = 0, size = 8)
  )

# Tables A5-A9: Models with alternative lag structures ----------------------------------------------

# Models without lags for FLFP
m0 <- lmer(log(F_M) ~ flfp + divorce_rate_lag1 + per_chr + (1 | party) + (1 | country_year), 
           data = gap_left, REML = FALSE)

m2 <- lmer(log(F_M) ~ flfp + gender_labor_equality + divorce_rate_lag1 + per_chr + share_women_lag1 + female_leader_lag1 + 
             (1 | party) + (1 | country_year), 
           data = gap_left, REML = FALSE)

m3 <- lmer(log(F_M) ~ flfp * gender_labor_equality + divorce_rate_lag1 + per_chr + share_women_lag1 + female_leader_lag1 + 
             (1 | party) + (1 | country_year), 
           data = gap_left, REML = FALSE)

order <- c("flfp", "gender_labor_equality", "share_women_lag1", "female_leader_lag1", "divorce_rate_lag1", "per_chr")

# Random effects for models without lags for FLFP
RE0 <- round(as.data.frame(VarCorr(m0))[, 4], 3)
RE2 <- round(as.data.frame(VarCorr(m2))[, 4], 3)
RE3 <- round(as.data.frame(VarCorr(m3))[, 4], 3)

stargazer(m0, m2, m3,
          title = "Models without lags for FLFP",
          add.lines = list(
            c("Var country-year", RE0[1], RE2[1], RE3[1]),
            c("Var party", RE0[2], RE2[2], RE3[2]),
            c("Var residual", RE0[3], RE2[3], RE3[3]),
            c("N country-years", nrow(ranef(m0)[[1]]), nrow(ranef(m2)[[1]]), nrow(ranef(m3)[[1]])),
            c("N parties", nrow(ranef(m0)[[2]]), nrow(ranef(m2)[[2]]), nrow(ranef(m3)[[2]]))
          ),
          order = order,
          star.cutoffs = c(0.05, 0.01, 0.001),
          dep.var.labels = 'log(F/M) party voter ratio',
          single.row = TRUE,
          type = "text")

# Models without lags for FLFP and divorce rate
m0 <- lmer(log(F_M) ~ flfp + divorce_rate + per_chr + (1 | party) + (1 | country_year), 
           data = gap_left, REML = FALSE)

m2 <- lmer(log(F_M) ~ flfp + gender_labor_equality + divorce_rate + per_chr + share_women_lag1 + female_leader_lag1 + 
             (1 | party) + (1 | country_year), 
           data = gap_left, REML = FALSE)

m3 <- lmer(log(F_M) ~ flfp * gender_labor_equality + divorce_rate + per_chr + share_women_lag1 + female_leader_lag1 + 
             (1 | party) + (1 | country_year), 
           data = gap_left, REML = FALSE)

order <- c("flfp", "gender_labor_equality", "share_women_lag1", "female_leader_lag1", "divorce_rate", "per_chr")

# Random effects for models without lags for FLFP and divorce rate
RE0 <- round(as.data.frame(VarCorr(m0))[, 4], 3)
RE2 <- round(as.data.frame(VarCorr(m2))[, 4], 3)
RE3 <- round(as.data.frame(VarCorr(m3))[, 4], 3)

stargazer(m0, m2, m3,
          title = "Models without lags for FLFP and divorce rate",
          add.lines = list(
            c("Var country-year", RE0[1], RE2[1], RE3[1]),
            c("Var party", RE0[2], RE2[2], RE3[2]),
            c("Var residual", RE0[3], RE2[3], RE3[3]),
            c("N country-years", nrow(ranef(m0)[[1]]), nrow(ranef(m2)[[1]]), nrow(ranef(m3)[[1]])),
            c("N parties", nrow(ranef(m0)[[2]]), nrow(ranef(m2)[[2]]), nrow(ranef(m3)[[2]]))
          ),
          order = order,
          star.cutoffs = c(0.05, 0.01, 0.001),
          dep.var.labels = 'log(F/M) party voter ratio',
          single.row = TRUE,
          type = "text")

# Models with two lags for FLFP
m0 <- lmer(log(F_M) ~ flfp_lag2 + divorce_rate_lag1 + per_chr + (1 | party) + (1 | country_year), 
           data = gap_left, REML = FALSE)

m2 <- lmer(log(F_M) ~ flfp_lag2 + gender_labor_equality + divorce_rate_lag1 + per_chr + share_women_lag1 + female_leader_lag1 + 
             (1 | party) + (1 | country_year), 
           data = gap_left, REML = FALSE)

m3 <- lmer(log(F_M) ~ flfp_lag2 * gender_labor_equality + divorce_rate_lag1 + per_chr + share_women_lag1 + female_leader_lag1 + 
             (1 | party) + (1 | country_year), 
           data = gap_left, REML = FALSE)

order <- c("flfp_lag2", "gender_labor_equality", "share_women_lag1", "female_leader_lag1", "divorce_rate_lag1", "per_chr")

# Random effects for models with two lags for FLFP
RE0 <- round(as.data.frame(VarCorr(m0))[, 4], 3)
RE2 <- round(as.data.frame(VarCorr(m2))[, 4], 3)
RE3 <- round(as.data.frame(VarCorr(m3))[, 4], 3)

stargazer(m0, m2, m3,
          title = "Models with two lags for FLFP",
          add.lines = list(
            c("Var country-year", RE0[1], RE2[1], RE3[1]),
            c("Var party", RE0[2], RE2[2], RE3[2]),
            c("Var residual", RE0[3], RE2[3], RE3[3]),
            c("N country-years", nrow(ranef(m0)[[1]]), nrow(ranef(m2)[[1]]), nrow(ranef(m3)[[1]])),
            c("N parties", nrow(ranef(m0)[[2]]), nrow(ranef(m2)[[2]]), nrow(ranef(m3)[[2]]))
          ),
          order = order,
          star.cutoffs = c(0.05, 0.01, 0.001),
          dep.var.labels = 'log(F/M) party voter ratio',
          single.row = TRUE,
          type = "text")

# Models with two lags for FLFP and divorce rates
m0 <- lmer(log(F_M) ~ flfp_lag2 + divorce_rate_lag2 + per_chr + (1 | party) + (1 | country_year), 
           data = gap_left, REML = FALSE)

m2 <- lmer(log(F_M) ~ flfp_lag2 + gender_labor_equality + divorce_rate_lag2 + per_chr + share_women_lag1 + female_leader_lag1 + 
             (1 | party) + (1 | country_year), 
           data = gap_left, REML = FALSE)

m3 <- lmer(log(F_M) ~ flfp_lag2 * gender_labor_equality + divorce_rate_lag2 + per_chr + share_women_lag1 + female_leader_lag1 + 
             (1 | party) + (1 | country_year), 
           data = gap_left, REML = FALSE)

order <- c("flfp_lag2", "gender_labor_equality", "share_women_lag1", "female_leader_lag1", "divorce_rate_lag2", "per_chr")

# Random effects for models with two lags for FLFP and divorce rates
RE0 <- round(as.data.frame(VarCorr(m0))[, 4], 3)
RE2 <- round(as.data.frame(VarCorr(m2))[, 4], 3)
RE3 <- round(as.data.frame(VarCorr(m3))[, 4], 3)

stargazer(m0, m2, m3,
          title = "Models with two lags for FLFP and divorce rate",
          add.lines = list(
            c("Var country-year", RE0[1], RE2[1], RE3[1]),
            c("Var party", RE0[2], RE2[2], RE3[2]),
            c("Var residual", RE0[3], RE2[3], RE3[3]),
            c("N country-years", nrow(ranef(m0)[[1]]), nrow(ranef(m2)[[1]]), nrow(ranef(m3)[[1]])),
            c("N parties", nrow(ranef(m0)[[2]]), nrow(ranef(m2)[[2]]), nrow(ranef(m3)[[2]]))
          ),
          order = order,
          star.cutoffs = c(0.05, 0.01, 0.001),
          dep.var.labels = 'log(F/M) party voter ratio',
          single.row = TRUE,
          type = "text")

# Models without lags for female leader
m0 <- lmer(log(F_M) ~ flfp_lag1 + divorce_rate_lag1 + per_chr + (1 | party) + (1 | country_year), 
           data = gap_left, REML = FALSE)

m2 <- lmer(log(F_M) ~ flfp_lag1 + gender_labor_equality + divorce_rate_lag1 + per_chr + share_women_lag1 + female_leader + 
             (1 | party) + (1 | country_year), 
           data = gap_left, REML = FALSE)

m3 <- lmer(log(F_M) ~ flfp_lag1 * gender_labor_equality + divorce_rate_lag1 + per_chr + share_women_lag1 + female_leader + 
             (1 | party) + (1 | country_year), 
           data = gap_left, REML = FALSE)

order <- c("flfp_lag1", "gender_labor_equality", "share_women_lag1", "female_leader", "divorce_rate_lag1", "per_chr")

# Random effects for models without lags for female leader
RE0 <- round(as.data.frame(VarCorr(m0))[, 4], 3)
RE2 <- round(as.data.frame(VarCorr(m2))[, 4], 3)
RE3 <- round(as.data.frame(VarCorr(m3))[, 4], 3)

stargazer(m0, m2, m3,
          title = "Models without lags for female leader",
          add.lines = list(
            c("Var country-year", RE0[1], RE2[1], RE3[1]),
            c("Var party", RE0[2], RE2[2], RE3[2]),
            c("Var residual", RE0[3], RE2[3], RE3[3]),
            c("N country-years", nrow(ranef(m0)[[1]]), nrow(ranef(m2)[[1]]), nrow(ranef(m3)[[1]])),
            c("N parties", nrow(ranef(m0)[[2]]), nrow(ranef(m2)[[2]]), nrow(ranef(m3)[[2]]))
          ),
          order = order,
          star.cutoffs = c(0.05, 0.01, 0.001),
          dep.var.labels = 'log(F/M) party voter ratio',
          single.row = TRUE,
          type = "text")





# Tables A10-A11: Models controlling for inflation ---------------------------------------------------

# Plot Figure: Consumer Price Index across time
ggplot(gap_left, aes(x = eyear, y = cpi, group = country, color = country)) +
  geom_line() +
  xlim(1980, 2031) +
  labs(title = "A. Consumer Price Index across time",
       x = "Year",
       y = "CPI") +
  theme_bw() +
  theme(legend.position = "none") +
  annotate("text", x = 1987, y = 13, label = "Portugal", color = "black", size = 3) +
  annotate("text", x = 1993, y = 17.5, label = "Greece", color = "black", size = 3) +
  annotate("text", x = 1985, y = 7, label = "Spain", color = "black", size = 3)

# Main model controlling for inflation (no lag)
m0 <- lmer(log(F_M) ~ flfp_lag1 + divorce_rate_lag1 + per_chr + cpi + 
             (1 | party) + (1 | country_year), 
           data = gap_left, REML = FALSE)

m2 <- lmer(log(F_M) ~ flfp_lag1 + gender_labor_equality + divorce_rate_lag1 + per_chr + cpi + 
             share_women_lag1 + female_leader_lag1 + 
             (1 | party) + (1 | country_year), 
           data = gap_left, REML = FALSE)

m3 <- lmer(log(F_M) ~ flfp_lag1 * gender_labor_equality + divorce_rate_lag1 + per_chr + cpi + 
             share_women_lag1 + female_leader_lag1 + 
             (1 | party) + (1 | country_year), 
           data = gap_left, REML = FALSE)

order <- c("flfp_lag1", "gender_labor_equality", "share_women_lag1", "female_leader_lag1", 
           "divorce_rate_lag1", "per_chr", "cpi")

# Extract and round random effects for each model (no lag)
RE0 <- round(as.data.frame(VarCorr(m0))[, 4], 3)
RE2 <- round(as.data.frame(VarCorr(m2))[, 4], 3)
RE3 <- round(as.data.frame(VarCorr(m3))[, 4], 3)

# Table (no lag)
stargazer(m0, m2, m3,
          title = "Models controlling for inflation",
          add.lines = list(
            c("Var country-year", RE0[1], RE2[1], RE3[1]),
            c("Var party", RE0[2], RE2[2], RE3[2]),
            c("Var residual", RE0[3], RE2[3], RE3[3]),
            c("N country-years", nrow(ranef(m0)[[1]]), nrow(ranef(m2)[[1]]), nrow(ranef(m3)[[1]])),
            c("N parties", nrow(ranef(m0)[[2]]), nrow(ranef(m2)[[2]]), nrow(ranef(m3)[[2]]))
          ),
          order = order,
          star.cutoffs = c(0.05, 0.01, 0.001),
          dep.var.labels = 'log(F/M) party voter ratio',
          single.row = TRUE,
          type = "text")

# Main model controlling for inflation (one year lag)
m0 <- lmer(log(F_M) ~ flfp_lag1 + divorce_rate_lag1 + per_chr + cpi_lag1 + 
             (1 | party) + (1 | country_year), 
           data = gap_left, REML = FALSE)

m2 <- lmer(log(F_M) ~ flfp_lag1 + gender_labor_equality + divorce_rate_lag1 + per_chr + cpi_lag1 + 
             share_women_lag1 + female_leader_lag1 + 
             (1 | party) + (1 | country_year), 
           data = gap_left, REML = FALSE)

m3 <- lmer(log(F_M) ~ flfp_lag1 * gender_labor_equality + divorce_rate_lag1 + per_chr + cpi_lag1 + 
             share_women_lag1 + female_leader_lag1 + 
             (1 | party) + (1 | country_year), 
           data = gap_left, REML = FALSE)

order <- c("flfp_lag1", "gender_labor_equality", "share_women_lag1", "female_leader_lag1", 
           "divorce_rate_lag1", "per_chr", "cpi_lag1")

# Extract and round random effects for each model (one year lag)
RE0 <- round(as.data.frame(VarCorr(m0))[, 4], 3)
RE2 <- round(as.data.frame(VarCorr(m2))[, 4], 3)
RE3 <- round(as.data.frame(VarCorr(m3))[, 4], 3)

# Table (one year lag)
stargazer(m0, m2, m3,
          title = "Models controlling for lagged inflation",
          add.lines = list(
            c("Var country-year", RE0[1], RE2[1], RE3[1]),
            c("Var party", RE0[2], RE2[2], RE3[2]),
            c("Var residual", RE0[3], RE2[3], RE3[3]),
            c("N country-years", nrow(ranef(m0)[[1]]), nrow(ranef(m2)[[1]]), nrow(ranef(m3)[[1]])),
            c("N parties", nrow(ranef(m0)[[2]]), nrow(ranef(m2)[[2]]), nrow(ranef(m3)[[2]]))
          ),
          order = order,
          star.cutoffs = c(0.05, 0.01, 0.001),
          dep.var.labels = 'log(F/M) party voter ratio',
          single.row = TRUE,
          type = "text")



# Tables A12-A13: Models controlling for social spending ---------------------------------------------------
# Plot Figure: Public Social Expenditure
ggplot(gap_left, aes(x = eyear, y = per_gdp_socexp, group = country, color = country)) +
  geom_line() +
  ylim(0, 40) +
  xlim(1980, 2027) +
  labs(title = "B. Public Social Expenditure",
       x = "Year",
       y = "% of GDP") +
  theme_bw() +
  theme(legend.position = "none") +
  annotate("text", x = 1983, y = 11, label = "Portugal", color = "black", size = 3)

# Main model controlling for %GDP on social expenditures
m0 <- lmer(log(F_M) ~ flfp_lag1 + divorce_rate_lag1 + per_chr + per_gdp_socexp + 
             (1 | party) + (1 | country_year), 
           data = gap_left, REML = FALSE)

m2 <- lmer(log(F_M) ~ flfp_lag1 + gender_labor_equality + divorce_rate_lag1 + per_chr + per_gdp_socexp + 
             share_women_lag1 + female_leader_lag1 + 
             (1 | party) + (1 | country_year), 
           data = gap_left, REML = FALSE)

m3 <- lmer(log(F_M) ~ flfp_lag1 * gender_labor_equality + divorce_rate_lag1 + per_chr + per_gdp_socexp + 
             share_women_lag1 + female_leader_lag1 + 
             (1 | party) + (1 | country_year), 
           data = gap_left, REML = FALSE)

order <- c("flfp_lag1", "gender_labor_equality", "share_women_lag1", "female_leader_lag1", 
           "divorce_rate_lag1", "per_chr", "per_gdp_socexp")

# Extract and round random effects for each model
RE0 <- round(as.data.frame(VarCorr(m0))[, 4], 3)
RE2 <- round(as.data.frame(VarCorr(m2))[, 4], 3)
RE3 <- round(as.data.frame(VarCorr(m3))[, 4], 3)

# Table
stargazer(m0, m2, m3,
          title = "Models controlling for social spending",
          add.lines = list(
            c("Var country-year", RE0[1], RE2[1], RE3[1]),
            c("Var party", RE0[2], RE2[2], RE3[2]),
            c("Var residual", RE0[3], RE2[3], RE3[3]),
            c("N country-years", nrow(ranef(m0)[[1]]), nrow(ranef(m2)[[1]]), nrow(ranef(m3)[[1]])),
            c("N parties", nrow(ranef(m0)[[2]]), nrow(ranef(m2)[[2]]), nrow(ranef(m3)[[2]]))
          ),
          order = order,
          star.cutoffs = c(0.05, 0.01, 0.001),
          dep.var.labels = 'log(F/M) party voter ratio',
          single.row = TRUE,
          type = "text")

# Main model controlling for %GDP on social expenditures last year
m0 <- lmer(log(F_M) ~ flfp_lag1 + divorce_rate_lag1 + per_chr + per_gdp_socexp_lag1 + 
             (1 | party) + (1 | country_year), 
           data = gap_left, REML = FALSE)

m2 <- lmer(log(F_M) ~ flfp_lag1 + gender_labor_equality + divorce_rate_lag1 + per_chr + per_gdp_socexp_lag1 + 
             share_women_lag1 + female_leader_lag1 + 
             (1 | party) + (1 | country_year), 
           data = gap_left, REML = FALSE)

m3 <- lmer(log(F_M) ~ flfp_lag1 * gender_labor_equality + divorce_rate_lag1 + per_chr + per_gdp_socexp_lag1 + 
             share_women_lag1 + female_leader_lag1 + 
             (1 | party) + (1 | country_year), 
           data = gap_left, REML = FALSE)

order <- c("flfp_lag1", "gender_labor_equality", "share_women_lag1", "female_leader_lag1", 
           "divorce_rate_lag1", "per_chr", "per_gdp_socexp_lag1")

# Extract and round random effects for each model
RE0 <- round(as.data.frame(VarCorr(m0))[, 4], 3)
RE2 <- round(as.data.frame(VarCorr(m2))[, 4], 3)
RE3 <- round(as.data.frame(VarCorr(m3))[, 4], 3)

# Table
stargazer(m0, m2, m3,
          title = "Models controlling for lagged social spending",
          add.lines = list(
            c("Var country-year", RE0[1], RE2[1], RE3[1]),
            c("Var party", RE0[2], RE2[2], RE3[2]),
            c("Var residual", RE0[3], RE2[3], RE3[3]),
            c("N country-years", nrow(ranef(m0)[[1]]), nrow(ranef(m2)[[1]]), nrow(ranef(m3)[[1]])),
            c("N parties", nrow(ranef(m0)[[2]]), nrow(ranef(m2)[[2]]), nrow(ranef(m3)[[2]]))
          ),
          order = order,
          star.cutoffs = c(0.05, 0.01, 0.001),
          dep.var.labels = 'log(F/M) party voter ratio',
          single.row = TRUE,
          type = "text")




# Tables A14-A15: Models controlling for SES ---------------------------------------------------

# Load the data
load("cses_respondent_party.RData")
cses_respondent_party <- dfidx(cses_respondent_party, idx = c("id", "parfam"))

# Baseline models
nw_women <- mlogit(voted_dummy ~ share_women_lag1 + female_leader_lag1 + gender_labor_equality |
                     education + age + hincome + divorced,
                   data = cses_respondent_party[cses_respondent_party$woman == 1 & cses_respondent_party$in_labor_force == 0, ],
                   reflevel = "30",
                   weights = dweight)  

w_women <- mlogit(voted_dummy ~ share_women_lag1 + female_leader_lag1 + gender_labor_equality |
                    education + age + hincome + divorced,
                  data = cses_respondent_party[cses_respondent_party$woman == 1 & cses_respondent_party$in_labor_force == 1, ],
                  reflevel = "30",
                  weights = dweight)  

nw_men <- mlogit(voted_dummy ~ share_women_lag1 + female_leader_lag1 + gender_labor_equality |
                   education + age + hincome + divorced,
                 data = cses_respondent_party[cses_respondent_party$woman == 0 & cses_respondent_party$in_labor_force == 0, ],
                 reflevel = "30",
                 weights = dweight) 

w_men <- mlogit(voted_dummy ~ share_women_lag1 + female_leader_lag1 + gender_labor_equality |
                  education + age + hincome + divorced,
                data = cses_respondent_party[cses_respondent_party$woman == 0 & cses_respondent_party$in_labor_force == 1, ],
                reflevel = "30",
                weights = dweight) 

stargazer(nw_women, w_women, nw_men, w_men,
          title = "Conditional Logistic Regressions of Vote Choice without Controlling for Religiosity",
          column.labels = c("Non-Working Women", "Working Women", "Non-Working Men", "Working Men"),
          star.cutoffs = c(0.05, 0.01, 0.001),
          dep.var.labels = 'Vote choice (Ref = Social Democratic party)',
          single.row = TRUE,
          type = "text")

# Models controlling for religiosity
nw_women <- mlogit(voted_dummy ~ share_women_lag1 + female_leader_lag1 + gender_labor_equality |
                     education + age + hincome + divorced + christian,
                   data = cses_respondent_party[cses_respondent_party$woman == 1 & cses_respondent_party$in_labor_force == 0, ],
                   reflevel = "30",
                   weights = dweight)  

w_women <- mlogit(voted_dummy ~ share_women_lag1 + female_leader_lag1 + gender_labor_equality |
                    education + age + hincome + divorced + christian,
                  data = cses_respondent_party[cses_respondent_party$woman == 1 & cses_respondent_party$in_labor_force == 1, ],
                  reflevel = "30",
                  weights = dweight)  

nw_men <- mlogit(voted_dummy ~ share_women_lag1 + female_leader_lag1 + gender_labor_equality |
                   education + age + hincome + divorced + christian,
                 data = cses_respondent_party[cses_respondent_party$woman == 0 & cses_respondent_party$in_labor_force == 0, ],
                 reflevel = "30",
                 weights = dweight) 

w_men <- mlogit(voted_dummy ~ share_women_lag1 + female_leader_lag1 + gender_labor_equality |
                  education + age + hincome + divorced + christian,
                data = cses_respondent_party[cses_respondent_party$woman == 0 & cses_respondent_party$in_labor_force == 1, ],
                reflevel = "30",
                weights = dweight) 

stargazer(nw_women, w_women, nw_men, w_men,
          title = "Conditional Logistic Regressions of Vote Choice",
          column.labels = c("Non-Working Women", "Working Women", "Non-Working Men", "Working Men"),
          star.cutoffs = c(0.05, 0.01, 0.001),
          dep.var.labels = 'Vote choice (Ref = Social Democratic party)',
          single.row = TRUE,
          type = "text")


# Table A16: Models controlling for incumbency status ------------------------------------------

m0 <- lmer(log(F_M) ~ flfp_lag1 + divorce_rate_lag1 + per_chr + incumbent_party + 
             (1 | party) + (1 | country_year), 
           data = gap_left, REML = FALSE)

m2 <- lmer(log(F_M) ~ flfp_lag1 + gender_labor_equality + divorce_rate_lag1 + per_chr + 
             share_women_lag1 + female_leader_lag1 + incumbent_party + 
             (1 | party) + (1 | country_year), 
           data = gap_left, REML = FALSE)

m3 <- lmer(log(F_M) ~ flfp_lag1 * gender_labor_equality + divorce_rate_lag1 + per_chr + 
             share_women_lag1 + female_leader_lag1 + incumbent_party + 
             (1 | party) + (1 | country_year), 
           data = gap_left, REML = FALSE)

order <- c("flfp_lag1", "gender_labor_equality", "share_women_lag1", "female_leader_lag1", 
           "divorce_rate_lag1", "per_chr", "incumbent_party")

# Extract and round random effects for each model
RE0 <- round(as.data.frame(VarCorr(m0))[, 4], 3)
RE2 <- round(as.data.frame(VarCorr(m2))[, 4], 3)
RE3 <- round(as.data.frame(VarCorr(m3))[, 4], 3)

# Table
stargazer(m0, m2, m3,
          title = "Models controlling for incumbency status",
          add.lines = list(
            c("Var country-year", RE0[1], RE2[1], RE3[1]),
            c("Var party", RE0[2], RE2[2], RE3[2]),
            c("Var residual", RE0[3], RE2[3], RE3[3]),
            c("N country-years", nrow(ranef(m0)[[1]]), nrow(ranef(m2)[[1]]), nrow(ranef(m3)[[1]])),
            c("N parties", nrow(ranef(m0)[[2]]), nrow(ranef(m2)[[2]]), nrow(ranef(m3)[[2]]))
          ),
          order = order,
          star.cutoffs = c(0.05, 0.01, 0.001),
          dep.var.labels = 'log(F/M) party voter ratio',
          single.row = TRUE,
          type = "text")



# Tables A17-A20: Models controlling for gender norms ---------------------------------------------------

# Load the data
load("ess_respondent.RData")

# Linear models
m1 <- lm(wmcpwrk ~ woman * pdwrk + age + income + education + divorced + christian + as.factor(country_essround), 
         data = ess_respondent, weights = anweight)

m2 <- lm(mnrgtjb ~ woman * pdwrk + age + income + education + divorced + christian + as.factor(country_essround), 
         data = ess_respondent, weights = anweight)

# Table
stargazer(m1, m2,
          title = "Work and Gender Norms",
          add.lines = list(c("Country-ESS round FE", "Yes", "Yes")),
          star.cutoffs = c(0.05, 0.01, 0.001),
          omit.stat = c("LL", "ser", "f"),
          single.row = TRUE,
          omit = "country_essround", 
          type = "text")




# Load the data
load("ess_respondent_party.RData")
ess_respondent_party <- dfidx(ess_respondent_party, idx = c("id", "parfam"))

# Baseline models
nw_women <- mlogit(voted_dummy ~ share_women_lag1 + female_leader_lag1 + gender_labor_equality |
                     education + age + income + divorced + christian,
                   data = ess_respondent_party[ess_respondent_party$woman == 1 & ess_respondent_party$pdwrk == 0, ],
                   reflevel = "30",
                   weights = anweight)  

w_women <- mlogit(voted_dummy ~ share_women_lag1 + female_leader_lag1 + gender_labor_equality |
                    education + age + income + divorced + christian,
                  data = ess_respondent_party[ess_respondent_party$woman == 1 & ess_respondent_party$pdwrk == 1, ],
                  reflevel = "30",
                  weights = anweight)  

nw_men <- mlogit(voted_dummy ~ share_women_lag1 + female_leader_lag1 + gender_labor_equality |
                   education + age + income + divorced + christian,
                 data = ess_respondent_party[ess_respondent_party$woman == 0 & ess_respondent_party$pdwrk == 0, ],
                 reflevel = "30",
                 weights = anweight) 

w_men <- mlogit(voted_dummy ~ share_women_lag1 + female_leader_lag1 + gender_labor_equality |
                  education + age + income + divorced + christian,
                data = ess_respondent_party[ess_respondent_party$woman == 0 & ess_respondent_party$pdwrk == 1, ],
                reflevel = "30",
                weights = anweight) 

stargazer(nw_women, w_women, nw_men, w_men,
          title = "Conditional Logistic Regressions of Vote Choice",
          column.labels = c("Non-Working Women", "Working Women", "Non-Working Men", "Working Men"),
          star.cutoffs = c(0.05, 0.01, 0.001),
          dep.var.labels = 'Vote choice (Ref = Social Democratic party)',
          single.row = TRUE,
          type = "text")



# Controlling for gender values: cut work for family
nw_women <- mlogit(voted_dummy ~ share_women_lag1 + female_leader_lag1 + gender_labor_equality |
                     wmcpwrk + education + age + income + divorced + christian,
                   data = ess_respondent_party[ess_respondent_party$woman == 1 & ess_respondent_party$pdwrk == 0, ],
                   reflevel = "30",
                   weights = anweight)  

w_women <- mlogit(voted_dummy ~ share_women_lag1 + female_leader_lag1 + gender_labor_equality |
                    wmcpwrk + education + age + income + divorced + christian,
                  data = ess_respondent_party[ess_respondent_party$woman == 1 & ess_respondent_party$pdwrk == 1, ],
                  reflevel = "30",
                  weights = anweight)  

nw_men <- mlogit(voted_dummy ~ share_women_lag1 + female_leader_lag1 + gender_labor_equality |
                   wmcpwrk + education + age + income + divorced + christian,
                 data = ess_respondent_party[ess_respondent_party$woman == 0 & ess_respondent_party$pdwrk == 0, ],
                 reflevel = "30",
                 weights = anweight) 

w_men <- mlogit(voted_dummy ~ share_women_lag1 + female_leader_lag1 + gender_labor_equality |
                  wmcpwrk + education + age + income + divorced + christian,
                data = ess_respondent_party[ess_respondent_party$woman == 0 & ess_respondent_party$pdwrk == 1, ],
                reflevel = "30",
                weights = anweight) 

stargazer(nw_women, w_women, nw_men, w_men,
          title = "Models Controlling for Gender Values (ESS waves 4 and 5). Higher Values of 'Disagreement to Cut Down Job' Mean More Egalitarian Gender Norms",
          column.labels = c("Non-Working Women", "Working Women", "Non-Working Men", "Working Men"),
          star.cutoffs = c(0.05, 0.01, 0.001),
          dep.var.labels = 'Vote choice (Ref = Social Democratic party)',
          single.row = TRUE,
          type = "text")


# Controlling for gender values: right to a job
nw_women <- mlogit(voted_dummy ~ share_women_lag1 + female_leader_lag1 + gender_labor_equality |
                     mnrgtjb + education + age + income + divorced + christian,
                   data = ess_respondent_party[ess_respondent_party$woman == 1 & ess_respondent_party$pdwrk == 0, ],
                   reflevel = "30",
                   weights = anweight)  

w_women <- mlogit(voted_dummy ~ share_women_lag1 + female_leader_lag1 + gender_labor_equality |
                    mnrgtjb + education + age + income + divorced + christian,
                  data = ess_respondent_party[ess_respondent_party$woman == 1 & ess_respondent_party$pdwrk == 1, ],
                  reflevel = "30",
                  weights = anweight)  

nw_men <- mlogit(voted_dummy ~ share_women_lag1 + female_leader_lag1 + gender_labor_equality |
                   mnrgtjb + education + age + income + divorced + christian,
                 data = ess_respondent_party[ess_respondent_party$woman == 0 & ess_respondent_party$pdwrk == 0, ],
                 reflevel = "30",
                 weights = anweight) 

w_men <- mlogit(voted_dummy ~ share_women_lag1 + female_leader_lag1 + gender_labor_equality |
                  mnrgtjb + education + age + income + divorced + christian,
                data = ess_respondent_party[ess_respondent_party$woman == 0 & ess_respondent_party$pdwrk == 1, ],
                reflevel = "30",
                weights = anweight) 

stargazer(nw_women, w_women, nw_men, w_men,
          title = "Models Controlling for Gender Values (ESS waves 4 and 5). Higher Values of 'Disagreement with Men Right to Jobs' Mean More Egalitarian Gender Norms",
          column.labels = c("Non-Working Women", "Working Women", "Non-Working Men", "Working Men"),
          star.cutoffs = c(0.05, 0.01, 0.001),
          dep.var.labels = 'Vote choice (Ref = Social Democratic party)',
          single.row = TRUE,
          type = "text")



# Tables A21-A22: Subgroup analyses by occupation groups -------------------------------


# Load the data
load("cses_respondent_party.RData")
cses_respondent_party <- dfidx(cses_respondent_party, idx = c("id", "parfam"))

# Women's class analysis
class <- subset(cses_respondent_party, cses_respondent_party$in_labor_force == 1)

prod <- mlogit(voted_dummy ~ share_women_lag1 + female_leader_lag1 + gender_labor_equality |
                 education + age + hincome + divorced + christian,
               data = class[class$woman == 1 & class$worker_type == "Production Workers", ],
               reflevel = "30",
               weights = dweight)  

service <- mlogit(voted_dummy ~ share_women_lag1 + female_leader_lag1 + gender_labor_equality |
                    education + age + hincome + divorced + christian,
                  data = class[class$woman == 1 & class$worker_type == "Service Workers", ],
                  reflevel = "30",
                  weights = dweight)  

clerk <- mlogit(voted_dummy ~ share_women_lag1 + female_leader_lag1 + gender_labor_equality |
                  education + age + hincome + divorced + christian,
                data = class[class$woman == 1 & class$worker_type == "Clerk Workers", ],
                reflevel = "30",
                weights = dweight)  

technical <- mlogit(voted_dummy ~ share_women_lag1 + female_leader_lag1 + gender_labor_equality |
                      education + age + hincome + divorced + christian,
                    data = class[class$woman == 1 & class$worker_type == "Technical Workers", ],
                    reflevel = "30",
                    weights = dweight)  

soccult <- mlogit(voted_dummy ~ share_women_lag1 + female_leader_lag1 + gender_labor_equality |
                    education + age + hincome + divorced + christian,
                  data = class[class$woman == 1 & class$worker_type == "Soccult Workers", ],
                  reflevel = "30",
                  weights = dweight)  

manager <- mlogit(voted_dummy ~ share_women_lag1 + female_leader_lag1 + gender_labor_equality |
                    education + age + hincome + divorced + christian,
                  data = class[class$woman == 1 & class$worker_type == "Manager Workers", ],
                  reflevel = "30",
                  weights = dweight)  

stargazer(prod, service, clerk, technical, soccult, manager,
          title = "Conditional Logistic Regressions of Women's Vote Choice by Occupation Group",
          column.labels = c("Production", "Service", "Clerk", "Technical", "Socio-Cultural", "Manager"),
          star.cutoffs = c(0.05, 0.01, 0.001),
          single.row = TRUE,
          dep.var.labels = 'Vote choice (Ref = Social Democratic party)',
          type = "text")

# Men's class analysis
prod <- mlogit(voted_dummy ~ share_women_lag1 + female_leader_lag1 + gender_labor_equality |
                 education + age + hincome + divorced + christian,
               data = class[class$woman == 0 & class$worker_type == "Production Workers", ],
               reflevel = "30",
               weights = dweight)  

service <- mlogit(voted_dummy ~ share_women_lag1 + female_leader_lag1 + gender_labor_equality |
                    education + age + hincome + divorced + christian,
                  data = class[class$woman == 0 & class$worker_type == "Service Workers", ],
                  reflevel = "30",
                  weights = dweight)  

clerk <- mlogit(voted_dummy ~ share_women_lag1 + female_leader_lag1 + gender_labor_equality |
                  education + age + hincome + divorced + christian,
                data = class[class$woman == 0 & class$worker_type == "Clerk Workers", ],
                reflevel = "30",
                weights = dweight)  

technical <- mlogit(voted_dummy ~ share_women_lag1 + female_leader_lag1 + gender_labor_equality |
                      education + age + hincome + divorced + christian,
                    data = class[class$woman == 0 & class$worker_type == "Technical Workers", ],
                    reflevel = "30",
                    weights = dweight)  

soccult <- mlogit(voted_dummy ~ share_women_lag1 + female_leader_lag1 + gender_labor_equality |
                    education + age + hincome + divorced + christian,
                  data = class[class$woman == 0 & class$worker_type == "Soccult Workers", ],
                  reflevel = "30",
                  weights = dweight)  

manager <- mlogit(voted_dummy ~ share_women_lag1 + female_leader_lag1 + gender_labor_equality |
                    education + age + hincome + divorced + christian,
                  data = class[class$woman == 0 & class$worker_type == "Manager Workers", ],
                  reflevel = "30",
                  weights = dweight)  

stargazer(prod, service, clerk, technical, soccult, manager,
          title = "Conditional Logistic Regressions of Men's Vote Choice by Occupation Group",
          column.labels = c("Production", "Service", "Clerk", "Technical", "Socio-Cultural", "Manager"),
          star.cutoffs = c(0.05, 0.01, 0.001),
          single.row = TRUE,
          dep.var.labels = 'Vote choice (Ref = Social Democratic party)',
          type = "text")



